PROJECT 1:

One-dimensional diffusion model for $CO_2$ absorption in the ocean

TMA4320 - Introduksjon til vitenskapelige beregninger

Group 3154

17.02.2021

---

Introduction

The concentration of $CO_2$ in the atmosphere has steadily increased since the industrial revolution leading to several serious climate related problems, and consequently an interest in understanding earth's changing climate. Earth's surfaces has the ability to absorb and store $CO_2$ from the atmosphere, playing a crucial role in the carbon cycle on earth. For instance, the global oceans acts as a natural carbon sink that help to buffer the emissions from human activity. This effect has great influence on the climate and it is therefore meaningful to construct scientific models that describes how the ocean reacts to atmospheric changes.

Quantifying the carbon absorbed by the global oceans is a complex problem and to simplify calculations everything is assumed constant in the horizontal direction, which means conditions only depend on depth and time. This is known as the water column model. Furthermore the project model employs the one-dimensional diffusion equation to describe how the concentration of dissolved inorganic carbon (DIC) in the ocean changes. Initially looking at shallow waters and constant atmospheric $CO_2$ concentration and subsequently in deeper water with an annual increase in atmospheric $CO_2$ concentration of $2.3$ppm.

---

Table of contents

1. Implementation
2. Libraries, constants and matrices

2.1 Libraries
2.2 Constants
2.3 Constructor for L, R and S matrices

3. Task 1 - Response to changing $CO_2$ concentration in shallow areas

3.1 Description
3.2 Functions
3.3 Results
3.4 Discussion

4. Task 2 - $CO_2$ absorption by the deep ocean

4.1 Description
4.2 Functions
4.3 Results
4.4 Discussion

5. Conclusion
6. References

---

1. Implementation

The first sections 2.1 Libraries and 2.2 Constants imports libraries, sets global constants. Values for depth, time and discretization steps for task 1 and task 2 are given in the first part of 3.2 Functions and 4.2 Functions.

Section 2.3 Constructor for L, R and S matrices creates the matrices $L$, $R$ and $S$ as described respectively in equations $(26a)$, $(26b)$ and $(24)$. The function is generalized to work for both task 1 and 2.

The computation of the DIC data for task 1 and 2 is in the sections 3.2 Functions and 4.2 Functions under the headline DIC concentration data. This task is split in two functions where one initiates the matrices and discretization, and calls the other function for the numerical computation. The functions are constructed this way such that the latter function can use the numba library (@njit) for increased efficiency. The time development of the DIC concentration is solved by the matrix equation $(28)$ which is a discretization of the one dimensional diffusion equation $(1)$. The time developed DIC concentration data is stored in a 2-dimensional array named C_data in the functions. Each element of C_data contains the concentration of DIC over depth and corresponds to a time value. Data for DIC, time and depth are stored in task_1 and task_2 for further visualization.

It is important to note that in task 1 and 2 the initial condition and development of $C_{eq}$ and consequently the $S$ matrix is different. In task 1 the zero concentration initial condition has been implemented as a concentration array filled with zeros, and the boundary conditions is given in the $S$ array with a constant $C_{eq}$ as in equation $(24)$. In task 2 the constant initial concentration is given as an array filled with the equilibrium value of $415$ppm. Since $C_{eq}$ varies with each time step in this problem, $S$ is now a matrix filled with the arrays from equation $(24)$ with varying $C_{eq}$.

Furthermore numerical integration of C_data is used to obtain the total mass of DIC and $CO_2$ in the ocean.

Note: $(number)$ represents the corresponding equation from project description [2].

---

2. Libraries, constants and matrices

2.1 Libraries

%matplotlib notebook

import numpy as np
from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy import linalg, sparse, integrate

from numba import njit

from IPython.display import set_matplotlib_formats, display, HTML, Image

# Plot parameters
fontsize = 12
plt.rcParams.update({
    "axes.titlesize":          fontsize + 5,
    "axes.labelsize":          fontsize,
    "axes.grid":               True,
    "axes.linewidth":          1,
    "lines.linewidth":         1.5,
    "lines.markersize":        7,
    "figure.figsize":          (9, 4.5),
    'axes.titlepad':           10,
    "ytick.labelsize":         fontsize,
    "xtick.labelsize":         fontsize,
    "ytick.major.pad":         5,
    "xtick.major.pad":         5,
    'legend.loc':              'best',
    "legend.fontsize":         fontsize,
    "legend.handlelength":     1.5,
    "legend.frameon":          True,
    "mathtext.fontset":        "stix",
    "font.family":             "STIXGeneral"
})
set_matplotlib_formats("svg")

2.2 Constants

# General
k_w            = 6.97e-5       # [m s^-1]            - Mass transfer coefficient
a_1            = 6.97e-7       # [s m^-1]            - Proportionality constant
u              = 10            # [m s^-1]            - Wind speed
H              = 5060          # [mol m^-3 atm^-1]   - Proportionality constant
P_CO2          = 415e-6        # [atm]               - Partial pressure CO2, assuming ideal gass
C_eq           = H * P_CO2     # [mol m^-3]          - DIC equilibrium concentration
startyear      = 2020          # [year]              - First year of simulations
ocean_area     = 360e12        # [m^2]               - Total area of ocean
C_atomic_mass  = 12            # [g mol^-1]          - Atomic mass of carbon


# Task 1
K_0_1          = 1e-3          # [m^2 s^-1]          - Constant 
K_a            = 2e-2          # [m^2 s^-1]          - Constant
K_b            = 5e-2          # [m^2 s^-1]          - Constant
z_a            = 7             # [m]                 - Constant
z_b            = 10            # [m]                 - Constant
L              = 100           # [m]                 - Depth of ocean


# Task 2
K_0_2          = 1e-4          # [m^2 s^-1]          - Constant
K_1_2          = 1e-2          # [m^2 s^-1]          - Constant
a_2            = 0.5           # [m^-1]              - Constant
z_0            = 100           # [m]                 - Constant

2.3 Constructor for L, R and S matrices

def matrix_constructor(
    depth,                            # (int)            - Total depth of ocean
    N,                                # (int)            - Discretization steps for depth
    delta_t,                          # (float)          - Timestep 
    matrix_name,                      # (string)         - Name of matrix to construct - L, R, S
    K_func,                           # (function)       - Function for diffusivity
    C_equilibrium = np.array([C_eq])  # (numpy.ndarray)  - Equilibrium concentration in ocean
):
    ### Discretization of diffusivity
    z_values, delta_z = np.linspace(0, depth, N + 1, retstep=True)
    K_values = K_func(z_values)
    
    
    ### Constants
    sign  = 1 if matrix_name == 'L' else -1 # L -> 1 | R -> -1
    alpha = delta_t / (2 * delta_z**2)
    gamma = 2 * alpha * k_w * delta_z * (1 - (K_values[1] - K_values[0]) / (2 * K_values[0]))
    
    
    ### S matrix
    if (matrix_name == 'S'):
        cont_column_vec = np.insert(np.zeros(N), 0, 2 * gamma)
        return np.dot(cont_column_vec[:, None], C_equilibrium[None, :]).T
    
    
    ### L and R matrix
    ## Main diagonal
    main_0 = 1 + sign * (2 * alpha * K_values[0] + gamma)
    main   = 1 + sign * 2 * alpha * K_values[1:]
    main   = np.insert(main, 0, main_0)
    
    ## Semidiagonals
    K_n_marked = K_values[2:] - K_values[:-2]
    
    # Upper diagonal
    upper_0 = -2 * alpha * K_values[0]
    upper   = -alpha / 4 * K_n_marked - alpha * K_values[1:-1]
    upper   = np.insert(upper, 0, upper_0) * sign
    
    # Lower diagonal
    lower_N = -2 * alpha * K_values[-1]
    lower   = alpha / 4 * K_n_marked - alpha * K_values[1:-1]
    lower   = np.append(lower, lower_N) * sign
    
    
    ### Total matrix
    return sparse.diags([upper, main, lower], offsets = [1, 0, -1]).toarray()

---

3. Task 1 - Response to changing $CO_2$ concentration in shallow areas

3.1 Description

Task 1 looks at the response in shallow areas of the ocean due to $CO_2$ changes in the atmosphere. The problem is of great importance because the increase in DIC concentration leads to acidification which is damaging for wildlife. The Norwegian continental shelf is such an area and the approximated depth is 100m.

The goal is to determine if the changes in the atmosphere and in the ocean are in sync, implying that the DIC concentration is in constant equilibrium with the atmosphere, or if it lags behind.

3.2 Functions

Initial values and data storage for task 1

task_1 = {
    'depth': {
        'value': 100,   # Depth of ocean
        'steps': 500    # Discretization steps
    },
    
    'days': {
        'value': 180,   # Days to simulate
        'steps': 800    # Discretization steps
    }
}

Diffusivity as a function of depth

def K_1(
    z # (float) - Depth from ocean top
):
    return K_0_1 + K_a * (z / z_a) * np.exp(-z / z_a) + K_b * (L - z) / z_b * np.exp(-(L - z) / z_b)

DIC concentration data

@njit
def get_C_data_1(
    z_N,          # (int)            - Depth discretization steps
    time_N,       # (int)            - Time discretization steps
    C,            # (numpy.ndarray)  - Initial DIC concentration
    L, R, S       # (numpy.ndarray)  - L, R, S matrices
):
    C_data = np.empty((time_N, z_N + 1))
    C_data[0, :] = C # First element
    
    for n in range(1, time_N):
        V = R.dot(C) + S
        C = np.linalg.solve(L, V)
        C_data[n, :] = C # Insert value
        
    return C_data

def task_1_DIC_data(
    depth,   # (int)  - Total depth of ocean
    z_N,     # (int)  - Discretization steps for depth
    days,    # (int)  - Time in days
    time_N   # (int)  - Discretization steps for time
):
    ### Discretization
    delta_t = days * 24 * 60 * 60 / (time_N - 1)
    time_values = np.linspace(0, days, time_N)
    z_values = np.linspace(0, depth, z_N + 1)
    
    ### Matrices
    L = matrix_constructor(depth, z_N, delta_t, 'L', K_1)
    R = matrix_constructor(depth, z_N, delta_t, 'R', K_1)
    S = matrix_constructor(depth, z_N, delta_t, 'S', K_1)[0]
    C = np.zeros(z_N + 1).T # Initial C
    
    ### Time development for C
    C_data = get_C_data_1(z_N, time_N, C, L, R, S)
    
    return C_data, time_values, z_values

Data visualization

def plot_DIC_C_time_dev_1(
    C_data,       # (np.ndarray)  - DIC concentration data
    time_values,  # (np.ndarray)  - Discretization of time
    time_N        # (int)         - Time discretization steps
):
    plt.figure()
    plt.title('Time development of DIC concentration')
    plt.ylabel('DIC concentration [$mol\cdot m^{-3}$]')
    plt.xlabel('Time [$days$]')
    
    # DIC equilibrium
    plt.plot(time_values, [C_eq] * time_N, label='$C_{eq}$')
    
    # DIC max
    plt.plot(time_values, C_data.max(axis=1), label='$C_{max}$')
    
    # DIC min
    plt.plot(time_values, C_data.min(axis=1), label='$C_{min}$')
    
    plt.legend()
    plt.show()

def plot_DIC_C_depth_anim(
    C_data,           # (np.ndarray)  - DIC concentration data
    z_values,         # (np.ndarray)  - Discretization of depth
    days,             # (int)         - Time in days
    depth,            # (int)         - Total depth of ocean
    anim_steps = 100, # (int)         - Animation steps
    tracing = False   # (boolean)     - Trace plot-line
):
    plt.ioff() # Turns off interactive plotting    
    
    fig, ax = plt.subplots()
    plt.title('DIC concentration as a function of depth')
    plt.xlabel('DIC concentration [$mol\cdot m^{-3}$]')
    plt.ylabel('Depth [$m$]')
    fig.set_tight_layout(True)
    ax.set_xlim(0, 2.5)  
    ax.set_ylim(0, 100)
    plt.gca().invert_yaxis()
    
    ## First line
    if tracing:
        lines = [ax.plot(C_data[0], z_values, color='b', alpha=.20, linewidth=1, label='Day: 0')[0]]
    else:
        lines = [ax.plot(C_data[0], z_values, color='b', linewidth=2, label='Day: 0')[0]]
    
    def update2(i):
        label = f'Day: {int(np.round(i / len(C_data) * days))}'
        
        if tracing: # See previous lines
            label0 = f'Day: 0 - {int(np.round(i / len(C_data) * days))}'
            
            ## Label for all previous lines
            lines[0].set_label(label0)
            
            ## Change previous line style
            lines[-1].set_label('_nolegend_')
            lines[-1].set_alpha(.20)
            lines[-1].set_color('b')
            lines[-1].set_linewidth(1)
        
            ## Add new line
            lines.append(ax.plot(C_data[int(i)], z_values, color='m', linewidth=2, label=label)[0])
        else: # Only one line
            lines[0].set_data(C_data[int(i)], z_values)
            lines[0].set_label(label)
        
        plt.legend(loc='upper right')
        
        return lines, ax
    
    ## Create animation
    anim = FuncAnimation(
        fig, update2, 
        frames=np.linspace(0, len(C_data) - 1, anim_steps), interval=100, 
        blit=True, repeat=False
    )
    
    ## Save animation
    return_anim = HTML(anim.to_jshtml())

    plt.close()
    plt.ion()
    
    return return_anim

Generate and store data

data = task_1_DIC_data(task_1['depth']['value'], task_1['depth']['steps'], task_1['days']['value'], task_1['days']['steps'])

task_1['C_data']      = data[0]
task_1['time_values'] = data[1]
task_1['z_values']    = data[2]

3.3 Results

Time development of DIC concentration

plot_DIC_C_time_dev_1(task_1['C_data'], task_1['time_values'], task_1['days']['steps'])

Plot shows extreme values of DIC concentration in the ocean as a function of time. The simulation shows development over 180 days, assumes zero initial concentration of DIC in the ocean and 415 ppm $CO_2$ in the atmosphere.

DIC concentration as a function of depth

display(plot_DIC_C_depth_anim(task_1['C_data'], task_1['z_values'], task_1['days']['value'], task_1['depth']['value'], 
                              tracing = True))

Once Loop Reflect

Animation shows how the concentration as a function of depth evolves through time.

3.4 Discussion

As illustrated in the animation above there is a gradual change in DIC concentration due to diffusivity. The change towards equilibrium happens at a faster rate in the shallow layer as a result of greater diffusivity. The animation also highlights that the concentration in the bottom part of the ocean shifts towards equilibrium at a much slower rate than the top. However, both plots still show that the concentration approaches equilibrium and reaches it within approximately 180 days. Considering that the concentration of $CO_2$ in the atmosphere varies only slowly from year to year, one could conclude that the changes are in fact in sync with the atmosphere, at least to a depth of 100 meters.

---

4. Task 2 - $CO_2$ absorption by the deep ocean

4.1 Description

The result of task 1 asserts that the DIC concentration in shallow water approaches atmospheric equilibrium quite fast. In light of this it is now appropriate to examine how the concentration of DIC changes in both shallow and deep water given a constant increase in the atmospheric $CO_2$ concentration over a decade. Additionally, the data generated during this examination will be used the see how the total mass of DIC in the ocean increases per year. This will then be compared to a more accurately calculated number, Gruber's number [1].

4.2 Functions

Initial values and data storage for task 2

task_2 = {
    'depth': {
        'value': 4000,  # Depth of ocean
        'steps': 4000   # Discretization steps
    },
    
    'years': {
        'value': 10,    # Years to simulate
        'steps': 120    # Discretization steps
    }
}

Diffusivity as a function of depth

def K_2(
    z # (float)  - Depth from ocean top
):
    return K_1_2 + (K_0_2 - K_1_2) / (1+ np.exp(-a_2 * (z - z_0)))

DIC equilibrium concentration

def C_eq_time(
    year # (float)  - years from 2020
):
    return H * (P_CO2 + 2.3e-6 * year)

DIC concentration data

@njit
def get_C_data_2(
    z_N,          # (int)            - Depth discretization steps
    time_N,       # (int)            - Time discretization steps
    C,            # (numpy.ndarray)  - Initial DIC concentration
    L, R, S       # (numpy.ndarray)  - L, R, S matrices
):
    C_data = np.empty((time_N, z_N + 1))
    C_data[0, :] = C # First element
    
    for i in range(1, time_N):
        V = R.dot(C) + .5 * (S[i - 1] + S[i])
        C = np.linalg.solve(L, V)
        C_data[i, :] = C # Insert value
        
    return C_data

def task_2_C_data(
    depth,   # (int)  - Total depth of ocean
    z_N,     # (int)  - Discretization steps for depth
    years,   # (int)  - Time in years
    time_N   # (int)  - Discretization steps for time
):
    ### Discretization
    delta_t = years * 365 * 24 * 60 * 60 / (time_N - 1)
    time_values = np.linspace(0, years, time_N)
    z_values, dz = np.linspace(0, depth, z_N + 1, retstep=True)
    
    ### Matrices
    L = matrix_constructor(depth, z_N, delta_t, 'L', K_2)
    R = matrix_constructor(depth, z_N, delta_t, 'R', K_2)
    S = matrix_constructor(depth, z_N, delta_t, 'S', K_2, C_eq_time(time_values))
    
    ### Initial C
    C = np.empty(z_N + 1)
    C.fill(C_eq)
    
    ### Data for C over time
    C_data = get_C_data_2(z_N, time_N, C, L, R, S)
    
    return C_data, time_values, z_values, dz

Carbon mass in ocean

def total_mass(
    C_data, # (np.ndarray)  - DIC concentration data
    dz      # (float)       - Discretization step size
):
    mass_data = np.array([])
    
    for C in C_data:
        # Total mass in ocean for given time
        mass = integrate.simps(C, dx=dz) * C_atomic_mass * ocean_area
        
        # Save total mass
        mass_data = np.append(mass_data, mass)
        
    return mass_data

$CO_2$ mass in ocean

def CO2_absorption(
    C_data, # (np.ndarray)  - DIC concentration data
    dz      # (float)       - Discretization step size
):
    # Difference in grams of carbon
    g_diff = (integrate.simps(C_data[-1], dx=dz) - integrate.simps(C_data[0], dx=dz)) * C_atomic_mass
    
    # Average over a decade for entire ocean
    CO2_diff_total = g_diff / task_2['years']['value'] * ocean_area
    
    return np.abs(CO2_diff_total)

Data visualization

def plot_DIC_C_time_dev_2_anim(
    C_data,            # (np.ndarray)  - DIC concentration data
    z_values,          # (np.ndarray)  - Discretization of depth
    years,             # (int)         - Time in days
    anim_steps = 100,  # (int)         - Animation steps
    tracing = False    # (boolean)     - Trace plot-line
):
    plt.ioff() # Turns off interactive plotting    
    
    fig, ax = plt.subplots()
    plt.title('DIC concentration as a function of depth')
    plt.xlabel('DIC concentration [$mol\cdot m^{-3}$]')
    plt.ylabel('Depth [$m$]')
    fig.set_tight_layout(True)
    ax.set_xlim(2.08, 2.22)
    plt.gca().invert_yaxis()
    
    # Initial graph
    line, = ax.plot(C_data[0], z_values, color='b', linewidth=2, label='Day: 0')
    
    def update2(i):
        label = f'Year: {int(startyear + np.round(i / len(C_data) * years))}'
        
        # Set new data and label
        line.set_data(C_data[int(i)], z_values)
        line.set_label(label)
        
        plt.legend(loc='lower right')
        
        return line, ax
     
    ## Create animation
    anim = FuncAnimation(
        fig, update2, 
        frames=np.linspace(0, len(C_data) - 1, anim_steps), interval=100, 
        blit=True, repeat=False
    )
    
    ## Save animation
    return_anim = HTML(anim.to_jshtml())

    plt.close()
    plt.ion()
    
    return return_anim

def plot_mass(
    C_data,      # (np.ndarray)  - DIC concentration data
    z_values,    # (np.ndarray)  - Discretization of depth
    dz,          # (float)       - Discretization step size
    time_values, # (np.ndarray)  - Discretization of time
    years        # (int)         - Time in years
):
    ### Get mass data
    mass_data = total_mass(C_data, dz)
    
    ### Plot data
    plt.figure()
    plt.title('Total DIC-mass in ocean')
    plt.xlabel('Year')
    plt.ylabel(r'DIC-mass [$10^{19} g$]')
    
    plt.plot(time_values + startyear, mass_data * 1e-19)
    
    plt.show()

Generate and store data

data = task_2_C_data(
    task_2['depth']['value'],
    task_2['depth']['steps'],
    task_2['years']['value'],
    task_2['years']['steps'] 
)

task_2['C_data']      = data[0]
task_2['time_values'] = data[1]
task_2['z_values']    = data[2] 
task_2['dz']          = data[3] 

4.3 Results

DIC concentration as a function of depth

display(plot_DIC_C_time_dev_2_anim(task_2['C_data'], task_2['z_values'], task_2['years']['value']))

Once Loop Reflect

Animation shows how the concentration as a function of depth evolves through time. Observe that the greatest change in concentration occurs in the upper layer of the ocean. This change reduces with increasing depth and in deep water there is hardly any change at all.

Total DIC-mass in ocean

plot_mass(task_2['C_data'], task_2['z_values'], task_2['dz'], task_2['time_values'], task_2['years']['value'])

The figure shows total mass of DIC in the entire ocean as a function of time. As expected, the mass of DIC in the water increases with an annual increase in atmospheric $CO_2$ concentration.

Total CO2 absorption in the entire ocean per year

ocean_CO2_absorption = CO2_absorption(task_2['C_data'], task_2['dz'])


print(
    'Carbon absorption in the entire global ocean per year in grams: ', 
    np.format_float_scientific(ocean_CO2_absorption, precision=3)
)

Carbon absorption in the entire global ocean per year in grams:  1.194e+16

4.4 Discussion

The animation of DIC concentration as a function of depth illustrates that in shallow water the concentration is synchronized with the atmospheric concentration due to turbulence. This is consistent with the results in task 1. Going deeper into the mixed layer the change flattens out but is still evident. In the mixed layer there is high diffusivity because the turbulence is largely involved in mixing the gas into the deeper water, as seen in the animation.

In the bottom part of the ocean the change in concentration is minimal, as expected with a diffusivity near zero.

The model yields $1.194 \cdot 10^{16}$g as the mass of DIC in the global ocean. Gruber's number [1], $2.5 \cdot 10^{15}$g, is in comparison 21% of the model's number. However the model used in this project is not flawless, for instance the oceans ability to absorb $CO_2$ is dependent on factors like location, temperature, ph and phytoplankton's absorption due too photosynthesis. All of these are ignored as a consequence of the water column model and the diffusion equation.

For instance, the disregard of the variations in the horizontal direction neglects that there might be local areas with net flow of $CO_2$ outwards. Another local variation which would obstruct the oceans ability to absorb $CO_2$ is that parts of the ocean is covered in ice, hindering surface absorption of $CO_2$. Taking this into account the model's number would decrease and approach Gruber's number.

---

5. Conclusion

The results from task 1 indicates that the DIC concentration in shallow waters is in sync with the $CO_2$ concentration in the atmosphere. The change towards equilibrium happens at a faster rate in the shallow layer in comparison to deeper levels, but both approach equilibrium within 180 days.

For task 2 this tendency persists in the shallow layer even with an annual increase of $CO_2$ in the atmosphere. In the mixed layer however, the DIC concentration increases but is not in sync with the atmospheric concentration. Furthermore in the deeper levels of the ocean the change in DIC concentration is negligible.

Based on our calculations the average absorption of $CO_2$ by the global oceans equals $1.194 \cdot 10^{16}$g, which is approximately 5 times larger than Gruber's number.

---

6. References

[1] N. Gruber, D. Clement, B. R. Carter, R. A. Feely, S. Van Heuven, M. Hoppema, M. Ishii, R. M. Key, A. Kozyr, S. K. Lauvset, et al. The oceanic sink for anthropogenic CO2 from 1994 to 2007. Science, 363(6432):1193–1199, 2019.

[2] T. Nordam. Project 1: One-dimensional diffusion model for CO2 absorption in the ocean. 2021.